Multispecies survey specifications

Model name	Full name	Latin name
Long_rough_dab	Long rough dab	Hippoglossoides platessoides
Green_halibut	Greenland halibut	Reinhardtius hippoglossoides
Mackerel	Mackerel	Scomber scombrus
Haddock	Haddock	Melongrammus aeglefinus
Saithe	Saithe	Pollachius virens
Redfish	Redfish	Sebastes mentella
Blue_whiting	Blue whiting	Micromesistius poutassou
Norwegian_ssh	Norwegian spring spawning herring	Clupea harengus
North_atl_cod	Northeast Atlantic cod	Gadus morhua
Polar_cod	Polar cod	Boreogadus saida
Capelin	Capelin	Mallotus villosus

Generating the census dataset (October 2020)

Note that survey config files must now name the survey.

Unchanging config files

NOBA2config.R looks like this (adjusted from Alfonso’s original):

d.name <- here("atlantisoutput","NOBA_march_2020")
functional.groups.file <- "nordic_groups_v04.csv"
biomass.pools.file <- "nordic_biol_v23.nc"
biol.prm.file <- "nordic_biol_incl_harv_v_007_3.prm"
box.file <- "Nordic02.bgm"
initial.conditions.file <- "nordic_biol_v23.nc"
run.prm.file <- "nordic_run_v01.xml"
scenario.name <- "nordic_runresults_01"
bioind.file <- "nordic_runresults_01BiomIndx.txt"
catch.file <- "nordic_runresults_01Catch.txt"
annage <- TRUE
fisheries.file <- "NoBAFisheries.csv"

omdimensions.R standardizes timesteps, etc. (this is part of atlantisom and should not need to be changed by the user):

#survey species inherited from omlist_ss
survspp <- omlist_ss$species_ss
# survey season and other time dimensioning parameters
# generalized timesteps all models
noutsteps <- omlist_ss$runpar$tstop/omlist_ss$runpar$outputstep
timeall <- c(0:noutsteps)
stepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutinc
midptyr <- round(median(seq(0,stepperyr)))

# model areas, subset in surveyconfig
allboxes <- c(0:(omlist_ss$boxpars$nbox - 1))

# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc

# survey selectivity (agecl based)
sp_age <- omlist_ss$funct.group_ss[, c("Name", "NumCohorts", "NumAgeClassSize")]

# should return all age classes fully sampled (Atlantis output is 10 age groups per spp)
n_age_classes <- sp_age$NumCohorts
# changed below for multiple species NOTE survspp alphabetical; NOT in order of fgs!!
# this gives correct names
age_classes <- sapply(n_age_classes, seq)
names(age_classes)<-sp_age$Name

n_annages <- sp_age$NumCohorts * sp_age$NumAgeClassSize
# changed below for multiple species
annages <- sapply(n_annages, seq)
names(annages)<-sp_age$Name

Change these survey and fishery config files

mssurvey_spring.R and mssurvey_fall.R configure the fishery independent surveys (in this census test, surveys sample all model polygons in all years and have efficiency of 1 for all species, with no size selectivity):

# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity

# Survey name
survey.name="BTS_spring_allbox_effic1"

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- 1 # spring survey

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

survtime <- survey_sample_full

# survey area
# should return all model areas
survboxes <- allboxes

# survey efficiency (q)
# should return a perfectly efficient survey 
surveffic <- data.frame(species=survspp,
                     efficiency=rep(1.0,length(survspp)))

# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(names(age_classes), each=n_age_classes),
#                     agecl=rep(c(1:n_age_classes),length(survspp)),
#                     selex=rep(1.0,length(survspp)*n_age_classes))

# for annage output uses names(annages) NOT alphabetical survspp
survselex <- data.frame(species=rep(names(annages), n_annages), #  
                        agecl=unlist(sapply(n_annages,seq)),
                        selex=rep(1.0,sum(n_annages)))

# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))

# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1

# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200

# Default survey configuration here has a range of efficiencies and selectivities
# To emulate a range of species in a single multispecies survey
# Also now happens in "spring" and "fall"
# Need to define survey season, area, efficiency, selectivity

# Survey name
survey.name="BTS_fall_allbox_effic1"

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- 3 # fall survey

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

survtime <- survey_sample_full

# survey area
# should return all model areas
survboxes <- allboxes

# survey efficiency (q)
# should return a perfectly efficient survey 
surveffic <- data.frame(species=survspp,
                     efficiency=rep(1.0,length(survspp)))

# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(survspp, each=n_age_classes),
#                     agecl=rep(c(1:n_age_classes),length(survspp)),
#                     selex=rep(1.0,length(survspp)*n_age_classes))

# for annage output
survselex <- data.frame(species=rep(names(annages), n_annages), #  
                        agecl=unlist(sapply(n_annages,seq)),
                        selex=rep(1.0,sum(n_annages)))

# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(100000, length(survspp)))

# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1

# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200

msfishery.R configures the fishery dependent data:

# Default fishery configuration here is a census
# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss 

#Number of years of data to pull
nyears <- 50

#Atlantis initialization period in years
burnin <- 30

# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc


# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample)  #total_sample defined in sardinesurvey.R
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[fstepperyr], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R

fishtime <- fish_times


# fishery sampling area
# should return all model areas, this assumes you see everything that it caught
fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))

# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))

#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr 

# fishery catch cv can be used in sample_survey_biomass
# perfect observation
fish_cv <- data.frame(species=survspp, cv=rep(0.01,length(survspp)))

These survey and fishery configurations were used previously to generate saved census survey output.

Read in previously generated survey and fishery outputs and plot them against true outputs

Using read_savedsurvs() we get back the same objects generated by the wrapper functions om_index() and om_comps().

source(here("config", "NOBA2config.R"))
omlist_ss <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
source(here("config/omdimensions.R"))

survObsBiom <- read_savedsurvs(d.name, 'survB')
age_comp_data <- read_savedsurvs(d.name, 'survAge')
len_comp_data <- read_savedsurvs(d.name, 'survLen')
wtage <- read_savedsurvs(d.name, 'survWtage')
annage_comp_data <- read_savedsurvs(d.name, 'survAnnAge')
annage_wtage <- read_savedsurvs(d.name, 'survAnnWtage')

# no function needed for saved fishery output as long as there is only one fishery specification file

catchbio_ss <- readRDS(file.path(d.name, paste0(scenario.name, "fishCatch.rds")))
catchbio_ss <- catchbio_ss[[1]] #in same time unit as surveys

fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]

#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- fish_age_comp$time/fstepperyr

fish_len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsLenComp.rds")))

#len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]

#add this to om_indices function so that this has years when read in
#fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))

fishwtage <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsWtAtAge.rds")))
fishwtage <- fishwtage[[1]] #need to change time units?

fish_annage_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_annage_comp <- fish_annage_comp[[1]]

#add this to om_indices function so that this has years when read in
#fish_annage_comp$time <- fish_annage_comp$time/fstepperyr

fish_annage_wtage <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullWtAtAge.rds")))
fish_annage_wtage <- fish_annage_wtage[[1]] #need to change time units?

Crunch a few things to get true age comps

# age classes
totNagecl <- omlist_ss$truenums_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

totN <- totNagecl %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueNagecl <- merge(totNagecl, totN)

# annual ages
totNage <- omlist_ss$truenumsage_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

#this should be exactly the same as totN but just making sure...
totN2 <- totNage %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueNage <- merge(totNage, totN2)

# fishery catch at age class
totCagecl <- omlist_ss$truecatchnum_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

totCAAcl <- totCagecl %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueCAAcl <- merge(totCagecl, totCAAcl)

# fishery catch at annual age
totCage <- omlist_ss$truecatchage_ss %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

totCAA <- totCage %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

trueCAA <- merge(totCage, totCAA)

My plotting functions: too specialized to include in the package?

# plot biomass time series facet wrapped by species
plotB <- function(dat, truedat=NULL){
  
    ggplot() +
    geom_line(data=dat, aes(x=time/stepperyr,y=atoutput, color="Survey Biomass"), 
              alpha = 10/10) +
    {if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True B"), alpha = 3/10)} + 
    theme_tufte() +
    theme(legend.position = "top") +
    xlab("model year") +
    ylab("tons") +
    labs(colour=scenario.name) +
    facet_wrap(~species, scales="free") 
  
}

# make a catch series function that can be split by fleet? this doesnt
# also note different time (days) from model timestep in all other output
plotC <- function(dat, truedat=NULL){
  
    ggplot() +
    geom_line(data=dat, aes(x=time/365,y=atoutput, color="Catch biomass"), 
              alpha = 10/10) +
    {if(!is.null(truedat)) geom_line(data=truedat, aes(x=time/365,y=atoutput, color="True Catch"), alpha = 3/10)} + 
    theme_tufte() +
    theme(legend.position = "top") +
    xlab("model year") +
    ylab("tons") +
    labs(colour=scenario.name) +
    facet_wrap(~species, scales="free") 
  
}

# plot length frequencies by timestep (one species)
plotlen <- function(dat, truedat=NULL){
  
    ggplot(dat, aes(upper.bins)) +
    geom_bar(aes(weight = atoutput)) +
    theme_tufte() +
    xlab("length (cm)") +
    ylab("number") +
    labs(subtitle = paste(scenario.name,
                          dat$species)) +
    facet_wrap(~time, ncol=6, scales="free_y")

}

# plot numbers at age by timestep (one species)
Natageplot <- function(dat, effN=1, truedat=NULL){
  ggplot() +
    geom_point(data=dat, aes(x=agecl, y=atoutput/effN, color="Est Comp")) +
    {if(!is.null(truedat)) geom_line(data=dat, aes(x=agecl, y=numAtAge/totN, color="True Comp"))} + 
    theme_tufte() +
    theme(legend.position = "bottom") +    
    xlab("age/agecl") +
    ylab("number") +
    labs(subtitle = paste(scenario.name,
                          dat$species)) + 
    facet_wrap(~time, ncol=6, scales="free_y")
}

# plot weight at age time series facet wrapped by species
wageplot <- function(dat, truedat=NULL){
  ggplot(dat, aes(time/stepperyr, atoutput)) +
    geom_line(aes(colour = factor(agecl))) +
    theme_tufte() +
    theme(legend.position = "bottom") +
    xlab("model year") +
    ylab("average individual weight (g)") +
    labs(subtitle = paste0(scenario.name)) +
    facet_wrap(c("species"), scales="free_y")
}

Visualize census surveys

We are using the full time series for biomass and weight at age, and a subsample of years 30-53 for compositions. In this test we see the effect of surveying in one season, with migratory species missed or only partially represented in some surveys, and with the overall surveyed biomass being at the lower or higher end of the range depending on the time of year. (Observation error is included in the survey biomass time series.) Samples for length and age assume a total of 100000 fish were randomly sampled for each species in each survey (surveffN defined in the survey config files).

Survey biomass

# compare with true output (all timesteps)
for(s in names(survObsBiom)){
  cat("  \n##### ",  s,"  \n")
  print(plotB(survObsBiom[[s]][[1]], omlist_ss$truetotbio_ss))
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

# plots survey only
# for(s in names(survObsBiom)){
#   print(plotB(survObsBiom[[s]][[1]]))
# }

Survey length composition

for(s in names(len_comp_data)){
  cat("  \n##### ",  s,"  \n")
  lcompsub <- as.data.frame(len_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
    group_by(species) %>%
    group_map(~ plotlen(.x), keep = TRUE)
  
  for(i in 1:length(lcompsub)) {
    print(lcompsub[[i]])
  }
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey age class composition (Atlantis age class)

for(s in names(age_comp_data)){
  cat("  \n##### ",  s,"  \n")
  acompsub <- as.data.frame(age_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
    group_by(species) %>%
    left_join(., trueNagecl) %>%
    #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
    group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp

  for(i in 1:length(acompsub)) {
    print(acompsub[[i]])
  }
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey age composition (annual ages)

for(s in names(annage_comp_data)){
  cat("  \n##### ",  s,"  \n")
  acompsub <- as.data.frame(annage_comp_data[[s]][[1]]) %>% filter(time %in% c(150:270)) %>%
    group_by(species) %>%
    left_join(., trueNage) %>%
    #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
    group_map(~ Natageplot(.x, effN = 100000, truedat = 1), keep = TRUE) # plots merged true age comp
  
  for(i in 1:length(acompsub)) {
    print(acompsub[[i]])
  }
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey weight at age (Atlantis age class)

for(s in names(wtage)){
  cat("  \n##### ",  s,"  \n")
  print(wageplot(wtage[[s]][[1]]))
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey iterpolated weight at age (annual ages)

# annage_wtage read in above 
# annage_wtage <- annage_wtage[[1]] #this still has second list component, diagnostic plot
# annage_wtage <- annage_wtage[[1]] #still in this old version but just removed from function.

for(s in names(annage_wtage)){
  cat("  \n##### ",  s,"  \n")
  print(wageplot(annage_wtage[[s]][[1]][[1]]))
  cat("  \n")
}

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Visualize perfectly observed fisheries

We are using the full time series for total catch and weight at age, and a subsample of years 30-35 for compositions. Fishery length and age composition is sampled every output timestep. Fishery effN for length and age sampling was 1000 fish per year (200 per time step), so this is reflected in the composition plots.

Fishery total catch

# compare with true catch, should match
plotC(catchbio_ss, omlist_ss$truecatchbio_ss)

Fishery length composition

lcompsub <- as.data.frame(fish_len_comp_data) %>% filter(time %in% c(150:175)) %>%
  group_by(species) %>%
  group_map(~ plotlen(.x), keep = TRUE)

for(i in 1:length(lcompsub)) {
  print(lcompsub[[i]])
}

Fishery catch at age (Atlantis age class)

acompsub <- as.data.frame(fish_age_comp) %>% filter(time %in% c(150:175)) %>%
  group_by(species) %>%
  left_join(., trueCAAcl) %>%
  #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
  group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp


for(i in 1:length(acompsub)) {
  print(acompsub[[i]])
}

Fishery catch at age (annual ages)

acompsub <- as.data.frame(fish_annage_comp) %>% filter(time %in% c(150:175)) %>%
  group_by(species) %>%
  left_join(., trueCAA) %>%
  #group_map(~ Natageplot(.x), keep = TRUE) # plots only sampled age comp
  group_map(~ Natageplot(.x, effN = 200, truedat = 1), keep = TRUE) # plots merged true age comp

for(i in 1:length(acompsub)) {
  print(acompsub[[i]])
}

Fishery weight at age (Atlantis age class)

wageplot(fishwtage)

Fishery interpolated weight at age (annual ages)

#diagnostic plot component still in this old version but just removed from function.
wageplot(fish_annage_wtage[[1]])

Multispecies survey specifications

Sarah Gaichas

28 October, 2020

Species surveyed

Generating the census dataset (October 2020)

Unchanging config files

Change these survey and fishery config files

Read in previously generated survey and fishery outputs and plot them against true outputs

Visualize census surveys

Survey biomass

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey length composition

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey age class composition (Atlantis age class)

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey age composition (annual ages)

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey weight at age (Atlantis age class)

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Survey iterpolated weight at age (annual ages)

BTS_fall_allbox_effic1

BTS_spring_allbox_effic1

Visualize perfectly observed fisheries

Fishery total catch

Fishery length composition

Fishery catch at age (Atlantis age class)

Fishery catch at age (annual ages)

Fishery weight at age (Atlantis age class)

Fishery interpolated weight at age (annual ages)